English

はじめに

このパッケージは多変量分析の結果を見やすく作図するためのパッケージです。 とりあえず、このパッケージはまだまだ開発中です。 バグは残っているでしょうし、足りない機能もたくさんあると思います。 一応注意はしていますが、もしかすると、研究をめちゃめちゃするような間違いがあるかもしれません。 今使えている機能が将来のバージョンで削除されることもあるかもしれませんし、使い方が変わるかもしれません。 と、いうことで、ご利用の際には注意してお使い下さい。 もしバグを見つけたり、こういう機能が欲しい、と思うことがあったりしたら、メール（アドレスは以下のコードをRにペースト：「rawToChar(as.raw(c(109, 111, 103, 64, 102, 102, 112, 114, 105, 46, 97, 102, 102, 114, 99, 46, 103, 111, 46, 106, 112)))」）もしくはGitHubにご連絡下さい。

インストール

以下のコマンドをRにコピー＆ペーストすれば 必要なパッケージ一式をインストールできるはずです。

install.packages(
    c("model.adapter", "partial.plot"), type = "source",
    repos = c(
        "http://florivory.net/R/repos", options()$repos
    )
)

クイックスタート（R魔法使い用）

まー、魔法使いの人は見ればわかるでしょう。 実例集も参考にしてください。

# データの読み込み。
data(iris)

# 予測モデルの作成。
model <- glm(
    Petal.Length ~ (Sepal.Length + Petal.Width) * Species, data = iris
)

# ライブラリの読み込み。
library(partial.plot)

# 萼片の長さと花弁の長さとの関係を作図。
info <- partial.plot(model, c("Sepal.Length", "Species"), pch = 16)

# 凡例を追加。
pp.legend(info, "topleft")

# 三次元プロット。
info <- partial.plot(
    model, c("Sepal.Length", "Petal.Width"), col = terrain.colors,
    theta = 20
)

基本的なつかいかた

データの準備とモデルの作成

今回はテスト用のデータにお約束のFisherのIrisデータを使います。 このデータは３種のアヤメ（setosa, versicolor, virginica）の萼片の長さ（Sepal.Length）、萼片の幅（Sepal.Width）、花弁の長さ（Petal.Length）、花弁の幅（Petal.Width）のデータが含まれています。

# データの読み込み。
data(iris)

# データの構造を見てみる。
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

さらに、種ごとに花弁の長さと萼片の長さ、萼片の幅、花弁の幅をプロットしてみます。

# データの可視化。
par(mfrow = c(2, 2))
plot(
    Petal.Length ~ Sepal.Length, data = iris,
    pch = as.numeric(Species) + 14, col = as.numeric(Species) + 1
)
plot(
    Petal.Length ~ Sepal.Width, data = iris,
    pch = as.numeric(Species) + 14, col = as.numeric(Species) + 1
)
plot(
    Petal.Length ~ Petal.Width, data = iris,
    pch = as.numeric(Species) + 14, col = as.numeric(Species) + 1
)

グラフを見ると、萼片の長さと花弁の長さの関係、萼片の幅と花弁の長さとの関係は種によって違いそうです。

このデータを用いて、花弁の長さを予測するモデルを作ることにします。 とりあえずモデルにはお約束の一般化線形混合モデル（GLM）を使い、説明変数には萼片の長さ、花弁の幅、種を用いることにします。 また、モデルに種と萼片の長さ、種と花弁の幅の交互作用を組み込んで、萼片の長さと花弁の長さ、花弁の幅と花弁の長さとの関係が種によって変わるかもしれない、ということを仮定します。

# 予測モデルの作成。
model <- glm(
    Petal.Length ~ (Sepal.Length + Petal.Width) * Species, data = iris
)

# 結果を確認します。
summary(model)
## 
## Call:
## glm(formula = Petal.Length ~ (Sepal.Length + Petal.Width) * Species, 
##     data = iris)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.72681  -0.11668   0.00073   0.12952   0.57611  
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     0.88128    0.47059   1.873  0.06318 .  
## Sepal.Length                    0.09342    0.09695   0.964  0.33693    
## Petal.Width                     0.45959    0.32429   1.417  0.15862    
## Speciesversicolor              -0.80183    0.60440  -1.327  0.18677    
## Speciesvirginica               -0.48261    0.60125  -0.803  0.42351    
## Sepal.Length:Speciesversicolor  0.32734    0.12315   2.658  0.00877 ** 
## Sepal.Length:Speciesvirginica   0.63569    0.11088   5.733 5.76e-08 ***
## Petal.Width:Speciesversicolor   0.80957    0.38007   2.130  0.03490 *  
## Petal.Width:Speciesvirginica   -0.28686    0.34738  -0.826  0.41033    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.05280405)
## 
##     Null deviance: 464.3254  on 149  degrees of freedom
## Residual deviance:   7.4454  on 141  degrees of freedom
## AIC: -4.7749
## 
## Number of Fisher Scoring iterations: 2

予測結果を見ると、Sepal.Length:Speciesversicolor、Sepal.Length:Speciesvirginica、Petal.Width:Speciesversicolorの推定値のP値が0.05以下なので、どうやら種によって萼片の長さと花弁の長さの関係、花弁の幅と花弁の長さの関係が異なっているようです。

モデルの予測結果の作図

上の結果を眺めていても、予測結果がどのような感じなのか、なかなかわかりません。 そこで、できあがったモデルを使ってモデルの予測結果を作図してみます。 基本的な使い方は以下の通りです。 以下の例ではモデルによって予測された種ごとの萼片の長さと花弁の長さとの関係を作図しています。

# ライブラリの読み込み。
library(partial.plot)

# 萼片の長さと花弁の長さとの関係を作図。
partial.plot(model, c("Sepal.Length", "Species"), pch = 16)

partial.plot()関数はまず一番目の引数にモデルの結果を入力します。 この場合はGLMの結果のオブジェクトmodelを使っています。

２番目の引数には結果を可視化する説明変数名を指定します。 この場合、萼片の長さを表す"Sepal.Length"と種を表す"Species"を指定しています。 因子型の変数は複数指定することができます。 以上の２つの引数は必ず指定する必要があります。

また、今回はグラフを見やすくするため、pch = 16を追加で指定し、 グラフのシンボルを変えています。

図に表された線はモデルによって予測された萼片の長さと花弁の長さの関係式、線の後ろの色の付いた部分は予測された関係式の信頼区間を、点は偏残差を表します。 偏残差は他の変数（この場合は花弁の長さ）で説明できなかった説明変数のばらつきを表します。

もし、花弁の幅と花弁の長さとの関係を作図したい場合、以下のように引数を指定します。

# 花弁の幅と花弁の長さとの関係を作図する。
partial.plot(model, c("Petal.Width", "Species"), pch = 16)

この場合、線は種ごとの花弁の幅と花弁の長さの関係式を、色の付いた部分は関係式の信頼区間を、点は偏残差を表します。

凡例の追加

ただ線や信頼区間、偏残差をプロットするだけでは色と種の対応関係がわかりません。 partial.plotの凡例追加機能を使って、凡例を追加することができます。

partial.plot()関数は図の設定などの情報を返します。例えば、

# 萼片の長さと花弁の長さとの関係を作図。
info <- partial.plot(model, c("Sepal.Length", "Species"), pch = 16)

というプログラムを書くと、図の情報がinfo変数に代入されます。 これを以下の例のようにpp.legend()関数に渡すと、できあがった図に凡例を追加することができます。

# 凡例を追加。
info <- partial.plot(model, c("Sepal.Length", "Species"), pch = 16)
pp.legend(info, "topleft")

pp.legend()関数の使い方はlegend()関数とほとんど同じです。 １番目の引数にはpartial.plot()の結果を指定します。 ２番目以降に入力した引数は全てlegend()関数にそのまま渡されます。

もう少し複雑なモデルの可視化

partial.plot()を使って、もう少し複雑なグラフを書くこともできます。 今度はRに含まれているCO2データセットを使います。 このデータセットはケベックとミシシッピからとってきた植物を使って低温処理実験を行い、その後、様々なCO2濃度でCO2吸収速度を測定したデータが入っています。

concが測定したCO2吸収速度、uptakeがCO2吸収速度、Treatmentが低温処理のありなし（あり：chilled、なし：nonchilled）、Typeが植物の由来（QuebecとMississippi）を表しています。

# CO2データを読み込む
data(CO2)

# データフレームに変換
class(CO2) <- "data.frame"

# データの傾向を見てみる。
plot(
    uptake ~ conc, data = CO2,
    col = c("blue3", "red3")[as.numeric(CO2$Type)],
    pch = c(16, 17)[as.numeric(CO2$Treatment)],
    ylim = c(0, 50)
)
legend(
    "bottomright",
    legend = c(
        "Mississippi - Chilled", "Mississippi - Non-Chilled",
        "Quebec - Chilled", "Quebec - Non-Chilled"
    ),
    col = rep(c("blue3", "red3"), each = 2),
    pch = c(16, 17, 16, 17)
)

この関係をglmでモデル化してみます。モデルの応答変数はCO2吸収速度（upteke）です。 上の図を見ると、CO2濃度とCO2吸収の関係は線形ではなさそうなので、説明変数にCO2濃度（conc）と二次の項（ I(conc ^ 2)）を加えることで二次式の当てはめを行います。 さらに、Type * Treatmentの項を加えることで、由来と低温処理の交互作用を見ることにします。

model <- glm(uptake ~ conc + I(conc^2) + Type * Treatment, data = CO2)
summary(model)
## 
## Call:
## glm(formula = uptake ~ conc + I(conc^2) + Type * Treatment, data = CO2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -11.8088   -3.2082    0.6907    3.0643    9.0688  
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       1.876e+01  1.829e+00  10.256 4.07e-16 ***
## conc                              6.666e-02  7.358e-03   9.060 8.27e-14 ***
## I(conc^2)                        -4.507e-05  6.579e-06  -6.851 1.51e-09 ***
## TypeMississippi                  -9.381e+00  1.472e+00  -6.373 1.20e-08 ***
## Treatmentchilled                 -3.581e+00  1.472e+00  -2.433  0.01728 *  
## TypeMississippi:Treatmentchilled -6.557e+00  2.082e+00  -3.150  0.00232 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 22.75242)
## 
##     Null deviance: 9707.0  on 83  degrees of freedom
## Residual deviance: 1774.7  on 78  degrees of freedom
## AIC: 508.63
## 
## Number of Fisher Scoring iterations: 2

どうやらモデル構築に使った変数全てがCO2吸収速度に影響していそうです。 ということで、例によって得られた関係式をpartial.plot()で可視化してみます。 このモデルで得られた関係性を可視化するには、以下のようにpartial.plot()を呼び出します。 ２番目の引数にc("conc", "Treatment", "Type")の３つを指定することでCO2濃度、処理、由来の３つの効果を同時に可視化することが可能です。

info <- partial.plot(
    model, c("conc", "Treatment", "Type"), data = CO2
)
pp.legend(info, "bottomright")

partial.plot()は複数の因子型によりグループごとの作図に対応しています。 ここではTreatmentとType、２つの因子型変数の組み合わせの４グループごとに作図を行いましたが、さらにたくさんの因子型引数でグループ分けすることも可能です。

三次元プロット

これまでの例ではモデルの予測結果を二次元に表していましたが、partial.plot()は三次元のプロットを描画することも可能です。 ここでも先ほどのFisherのIrisデータを使います。 今回は（テスト用に）種の違いを無視して花弁の長さを萼片の長さと花弁の幅で説明してみます。

# データの読み込み。
data(iris)

# 予測モデルの作成。
model <- glm(
    Petal.Length ~ Sepal.Length + Petal.Width, data = iris
)

# 予測結果を見る。
summary(model)
## 
## Call:
## glm(formula = Petal.Length ~ Sepal.Length + Petal.Width, data = iris)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.15506  -0.21920  -0.02115   0.25986   1.35204  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.50714    0.33696  -4.473 1.54e-05 ***
## Sepal.Length  0.54226    0.06934   7.820 9.41e-13 ***
## Petal.Width   1.74810    0.07533  23.205  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1625972)
## 
##     Null deviance: 464.325  on 149  degrees of freedom
## Residual deviance:  23.902  on 147  degrees of freedom
## AIC: 158.18
## 
## Number of Fisher Scoring iterations: 2

予測結果を見ると萼片の長さも花弁の花弁の長さに有意な影響を与えているようです。 そこで、partial.plot()を使って予測結果を可視化してみます。

partial.plot(
    model, c("Sepal.Length", "Petal.Width"), col = terrain.colors, theta = 20
)

デフォルトでは花弁の長さと萼片の長さ、花弁の幅の関係がpersp()関数を使って描画されます。 三次元プロットに使う関数を変えるにはpartial.plotのfun.3d引数に使いたい関数を指定します。 現在、persp(), image()、contour()での描画に対応しています。

partial.plot(
    model, c("Sepal.Length", "Petal.Width"), fun.3d = image,
    col = terrain.colors
)

partial.plot(
    model, c("Sepal.Length", "Petal.Width"), fun.3d = contour, col = "black"
)

もし、お使いのパソコンにrglパッケージがインストールしてあるなら、persp3d()関数も使うことができます。 ブラウザによっては結果が表示できないので、試してみたい場合にはrglパッケージをインストールして、お使いのRで試してみて下さい。

# rglパッケージをインストール。
install.packages("rgl")

# 三次元プロット
partial.plot(
    model, c("Sepal.Length", "Petal.Width"), fun.3d = persp3d,
    col = terrain.colors
)

グラフの見た目の変更

いくつかの設定を変えることで、グラフの見た目を変えることができます。

作図する要素の変更

draw.residual引数をFALSEにすると、残差のシンボルを描画しないことができます。

info <- partial.plot(
    model, c("Sepal.Length", "Species"), draw.residual = FALSE
)
pp.legend(info, "topleft")

また、draw.relationship引数をFALSEにすると、モデルによって予測された関係を描画しないことができます。

info <- partial.plot(
    model, c("Sepal.Length", "Species"), draw.relationship = FALSE,
    pch = 16
)
pp.legend(info, "topleft")

さらに、draw.interval引数をFALSEにすると、予測された関係の信頼区間を描画しないこともできます。

info <- partial.plot(
    model, c("Sepal.Length", "Species"), draw.interval = FALSE,
    pch = 16
)
pp.legend(info, "topleft")

extrapolate引数をTRUEにすると、説明変数の範囲外の値も含めた予測範囲を描画します。

info <- partial.plot(
    model, c("Sepal.Length", "Species"), pch = 16, extrapolate = TRUE
)
pp.legend(info, "topleft")

色の変更

rainbow()関数やheat.colors()関数など、一番目の引数に色の数を渡すと色を返す関数を使ってグラフの色を変えることができます。 ?rainbowで出てくるヘルプに載っている関数などが使えます。 個人的にはviridisパッケージのviridis()関数で精製される色がおすすめです。

# グラフの色を変える。
info <- partial.plot(
    model, c("Sepal.Length", "Species"), pch = 16, col = rainbow
)
# 凡例の色も自動調節される。
pp.legend(info, "topleft")

また、名前付きの色ベクトルを作ることで、個別に色を指定することもできます。 以下の例のようにc()関数を用いてグループ名 = 色という形式で色ベクトルを作り、それをcol引数に渡すことで各グループ（この場合は種）の色を指定することができます。

# 色ベクトルを用意。
col <- c(
    setosa = "darkgreen", versicolor = "blue", virginica = "orange2"
)
# 作図
info <- partial.plot(
    model, c("Sepal.Length", "Species"), pch = 16, col = col
)
pp.legend(info, "topleft")

プロットのシンボル、線の種類、線の太さの変更

色の指定と同じ方法で、プロットのシンボル、線の種類、線の太さを変更することができます。

シンボルの変更

info <- partial.plot(
    model, c("Sepal.Length", "Species"), pch = 16:18, col = rainbow
)
pp.legend(info, "topleft")

線の種類の変更

info <- partial.plot(
    model, c("Sepal.Length", "Species"), lty = c("solid", "dashed", "dotted"),
    col = "black"
)
pp.legend(info, "topleft")

線の太さの変更

info <- partial.plot(
    model, c("Sepal.Length", "Species"), lwd = c(1, 4, 8), col = "black"
)
pp.legend(info, "topleft")

グラフの上書き

add = TRUEを指定すると、既存のグラフの上にグラフを上書きすることができます。

model <- glm(
    Petal.Length ~ (Sepal.Length + Petal.Width) * Species, data = iris
)
plot(Petal.Length ~ Sepal.Length, data = iris, pch = 16)
partial.plot(model, c("Sepal.Length", "Species"), add = TRUE)

ラベル区切り文字の変更

複数の因子型説明変数を描画するとき、レジェンドで因子名の間を区切る文字を変えることができます。

model <- glm(uptake ~ conc + I(conc^2) + Type * Treatment, data = CO2)
info <- partial.plot(
    model, c("conc", "Treatment", "Type"), data = CO2, sep = "/"
)
pp.legend(info, "bottomright")

信頼区間の変更

線の周りに描画される信頼区間はデフォルトでは95%信頼区間になっています。 interval.levels引数に0から1の値を与えることで、この値も変更することが可能です。

par(mfrow = c(2, 2))
model <- glm(uptake ~ conc + I(conc^2) + Type * Treatment, data = CO2)
# 95%（デフォルト）
info <- partial.plot(
    model, c("conc", "Treatment", "Type"), data = CO2,
    main = "interval.levels = 0.95"
)
# 80%
info <- partial.plot(
    model, c("conc", "Treatment", "Type"), data = CO2, interval.levels = 0.8,
    main = "interval.levels = 0.8"
)
# 70%
info <- partial.plot(
    model, c("conc", "Treatment", "Type"), data = CO2, interval.levels = 0.7,
    main = "interval.levels = 0.7"
)
# 60%
info <- partial.plot(
    model, c("conc", "Treatment", "Type"), data = CO2, interval.levels = 0.6,
    main = "interval.levels = 0.6"
)

par(mfrow = c(1, 1))

リンク関数

GLMやGLMMなどのモデルを用いると、リンク関数が適用された応答変数に対して、線形モデルの当てはめが行われます。 GLM全般を説明するのは大変なので、詳細はGoogle先生に聞いてみて下さい。 例えば、下の例のように分布にガンマ分布、リンク関数にlogを用いると、GLMは \[ Petal.Length = exp(切片 + 係数１ \times Sepal.Length + 係数２ \times Petal.Width) \] という関係式を推定します。

model <- glm(
    Petal.Length ~ (Sepal.Length + Petal.Width), data = iris,
    family = Gamma(log)
)
summary(model)
## 
## Call:
## glm(formula = Petal.Length ~ (Sepal.Length + Petal.Width), family = Gamma(log), 
##     data = iris)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.39018  -0.13095  -0.01123   0.11138   0.40059  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.19569    0.15091  -1.297 0.196754    
## Sepal.Length  0.10529    0.03106   3.390 0.000897 ***
## Petal.Width   0.64334    0.03374  19.069  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.03261211)
## 
##     Null deviance: 44.6546  on 149  degrees of freedom
## Residual deviance:  4.8219  on 147  degrees of freedom
## AIC: 271.38
## 
## Number of Fisher Scoring iterations: 5

モデルの推定結果で切片に対応するのはsummary(model)の結果の(Intercept)のEstimate（-0.19569）、係数１に対応するのはSepal.LengthのEstimate（0.10529）、係数２に対応するのはPetal.WidthのEstimate（0.64334）です。

partial.plot()は作図の際、応答変数と同じ単位（普通の単位）・線形予測子の単位（リンク関数を適用したあとの単位）、両方の単位で作図を行うことができます。 作図の単位を制御するにはtype引数を用い、type引数に"response"を指定すると応答変数の単位で、type引数にlinkを指定すると線形予測子の単位で作図を行います。

partial.plot(
    model, c("Sepal.Length", "Petal.Width"), type = "response",
    col = terrain.colors, main = "response"
)

partial.plot(
    model, c("Sepal.Length", "Petal.Width"), type = "link",
    col = terrain.colors, main = "link"
)

識別モデルの作図

partial.plotを用いると、応答変数が連続値の回帰モデルだけではなく、識別モデルの予測確率も作図することができます。 識別モデルの予測値を計算するには、今のところtype = "prob"を指定する必要があります。 将来的には自動でモデルの種類を判定して作図するようにする予定です。 また、識別モデルでは残差の作図はサポートされていません。

library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
model <- randomForest(Species ~ ., data = iris)
partial.plot(
    model, "Petal.Length", positive.class = "setosa", type = "prob",
    col = "red", resolution = 20, draw.residual = FALSE
)
partial.plot(
    model, "Petal.Length", positive.class = "versicolor", add = TRUE,
    type = "prob", col = "blue", resolution = 20, draw.residual = FALSE
)
partial.plot(
    model, "Petal.Length", positive.class = "virginica", add = TRUE,
    type = "prob", col = "green", resolution = 20, draw.residual = FALSE
)

ヒストグラムの追加

試験的に説明変数のヒストグラムを描く機能を追加してあります。

model <- glm(Petal.Length ~ ., data = iris)
info <- partial.plot(
    model, c("Sepal.Length", "Species"), draw.hist = TRUE
)

その他の設定

partial.plot()関数は普通のplot()関数と同じように、グラフィックパラメーターを変えて、グラフの見た目を変えることができます。

# モデル作成
model <- glm(
    Petal.Length ~ (Sepal.Length + Petal.Width) * Species, data = iris
)
# シンボルのサイズ、軸ラベルも変えてみる。
info <- partial.plot(
    model, c("Sepal.Length", "Species"), cex = 1.5,
    xlab = "Sepal length (mm)", ylab = "Petal Length"
)
pp.legend(info, "topleft")

結果の再利用

モデルによっては予測値を計算するのに計算時間がかかることがあります。 partial.plot()は結果に計算済みの関係式や偏残差情報を記憶しているので描画の設定だけを変えて、グラフを書き直すことができます。

par(mfrow = c(2, 2))

# 最初のpartial.plotを実行。
info <- partial.plot(model, c("Sepal.Length", "Species"), pch = 16)
pp.legend(info, "topleft", cex = 0.5)

# 結果を再利用して、シンボルと関係式描画の設定を変える。
info <- partial.plot(info, pch = 6, draw.relationship = FALSE)
pp.legend(info, "topleft", cex = 0.5)

# 結果を再利用して、シンボル、色、信頼区間描画の設定を変える。
info <- partial.plot(info, col = rainbow, draw.interval = FALSE)
pp.legend(info, "topleft", cex = 0.5)

# 結果を再利用して、シンボル、色、偏残差描画の設定を変える。
info <- partial.plot(info, col = heat.colors, draw.residual = FALSE)
pp.legend(info, "topleft", cex = 0.5)

par(mfrow = c(1, 1))

現在、書き換えることのできる設定はfun.3d, draw.residual, draw.relationship, draw.interval, col, xlab, ylab, zlab, その他のグラフィックパラメーターです。

レジェンドの複数因子間の区切り文字sepは実装上の都合により変更することができません。 要望があれば、そのうち修正するかもしれません。

また、予測された関係式の信頼区間を指定するinterval.levelsも同様の理由により変更できません。 こちらはどうしても再計算が必要なので、修正が難しいと思います。

対応しているモデル

今のところ、lm()、glm()、glm.nb()、lme()、lmer()、glmer()、glmer.nb()、glmmadmb()、MCMCglmm()、cforest()、ctree()、svm()、randomForest()、ranger()、rpart()、tree()に対応しています（たぶん）。 各モデルへの対応はmodel.adapterパッケージに依存しているので、 将来的にはmodel.adapterパッケージが対応する他のモデルにも対応する予定です。

計算の詳細

partial.plotでは説明変数と応答変数の関係の推定にemmeansパッケージを使用しています。 もし対象のモデルにemmeansパッケージが対応している場合、説明変数と応答変数の関係はemmeansで計算されます。 今のところ、lm、glm、glm.nb、lme、lmer、glmer、glmer.nb、glmmadmb、MCMCglmmの計算にemmeansが使われます。 モデルにemmeansが対応している場合、対象の説明変数以外の変数は平均値に固定され、予測値が計算されます。

対象のモデルにemmeansが対応していない場合、計算はpartial.plotによって行われます。 この場合、対象の説明変数以外の説明変数の値には元のデータセット全体の値が使われ、予測値の平均値を用いて説明変数と応答変数の関係が描画されます。

モデルごとの注意点

glmer

2021-06-23現在、モデル作成時にデータをscale関数で標準化すると、データフレーム内の列がmatrixになるため、predict関数がエラーになり、partial.plotの実行が止まるようです。 このエラーを回避するため、データを標準化するときには

data$column1 <- c(scale(data$column1))

のような感じで標準化した結果をベクトルに変換してください。

glmer.nb

モデル式に関数（’offset()’や’log()’など）が含まれていると、データの取得がうまくいかないため、作図を行うことができないようです。 このようなモデル式を用いたglmer.nbで作図を行う時には、data引数に解析に使ったデータを指定してください。

# 例
model <- glmer.nb(y ~ log(x) + (1 | random), data = dat)
partial.plot(model, "x", data = dat)

offset項のあるモデル

offset項を用いた予測値の計算

partial.plotではoffset項は以下のように扱われます。

予測値と信頼区間
他の説明変数と同様、offset項に使われた変数も平均値が予測値の計算には使われます。例えば応答変数が個体数、offset項が調査面積の場合、予測値と信頼区間は調査面積の平均値を使って計算されます。
偏残差
offset項に指定された変数の値には、モデル作成に使われた値をそのまま用いて予測値を計算し、実際の応答変数の値との差分から偏残差を計算します。

このため、そのまま作図を行うと、偏残差と予測値が対応しないことがあります。対応方法は２つあり、１つはdraw.residual = FALSEオプションを用いて偏残差を描画しない、という方法です。もう１つは以下の例のようにあらかじめoffset項あたりの応答変数のデータを準備し、partial.plotを実行する方法です。

offset項非対応のモデル

現在、以下のモデルはオフセット項のある予測に対応していません。

  • glmmADMB::glmmadmb
#-----------------------------------------------------------------------------
#   データとモデルの準備
#-----------------------------------------------------------------------------

# MASSパッケージから車の保険請求データを準備。
utils::data(Insurance, package="MASS")

# 年齢をカテゴリーから数値データに変換。
Insurance$Age.numeric <- c(25, 27, 32.5, 35)[as.numeric(Insurance$Age)]

# 保険請求数を地区・車の種類・年齢でモデル化。
# offset項は契約者数。
model <- glm(
    Claims ~ District + Group + Age.numeric + offset(log(Holders)),
    data = Insurance, family = poisson
)

#-----------------------------------------------------------------------------
#   作図１：失敗例
#-----------------------------------------------------------------------------

# 平均の契約者数を使ってpartial.plotを描画。
# そのまま描画すると偏残差と予測値が対応しない。
pp <- partial.plot(model, c("Age.numeric", "Group"))
pp.legend(pp, "topright")

#-----------------------------------------------------------------------------
#   対策１：偏残差を描画しない
#   この場合は平均契約者数あたりの請求数が描画される。
#-----------------------------------------------------------------------------
pp <- partial.plot(model, c("Age.numeric", "Group"), draw.residual = FALSE)
pp.legend(pp, "topright")

#-----------------------------------------------------------------------------
#   対策２：offset項あたりの応答変数を用いる
#   この場合は契約者数あたりの請求数を先に計算し、
#   契約者数を１に固定したデータを用いる。
#-----------------------------------------------------------------------------

# 描画用にデータのコピーを作成。
Insurance2 <- Insurance

# 応答変数を契約者あたりの請求数に変換。
Insurance2$Claims <- Insurance2$Claims / Insurance2$Holders

# 契約者数を１に固定。
Insurance2$Holders <- 1

# 契約者数あたりの請求数に対してモデルの予測値を描画。
pp <- partial.plot(model, c("Age.numeric", "District"), data = Insurance2)
pp.legend(pp, "topright")

既知の問題

今のところ、以下の問題があることがわかっています。

おまけ

作図に使っている以下の関数も何かの役立つかもしれません。

color.ramp()

color.ramp()関数は因子型変数から色を表す文字列ベクトルを生成します。 以下の例のように使います。

library(partial.plot)
plot(
    Petal.Length ~ Sepal.Length, col = color.ramp(Species),
    data = iris, pch = 16
)

partial.plotと同じように色を変えることができます。

plot(
    Petal.Length ~ Sepal.Length, col = color.ramp(Species, rainbow),
    data = iris, pch = 16
)

col <- c(
    setosa = "darkgreen", versicolor = "blue", virginica = "orange2"
)
plot(
    Petal.Length ~ Sepal.Length, col = color.ramp(Species, col),
    data = iris, pch = 16
)

switch.par()

color.ramp()と同じようにグラフィックパラメーターを変更することができます。

plot(Petal.Length~Sepal.Length, data = iris, pch = switch.par(Species, 16:18))

gg.colors()

ggplot()と似たような色を生成します。

barplot(1:10, col = gg.colors(10))

更新履歴